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ABSTRACT 


There  are  basically  two  techniques  used  to  solve  the  Navier-Stokes 
equations  for  fluid  flow.  These  techniques  are  implicit  and  explicit 
methods.  Both  methods  along  with  characteristic  boundary  conditions 
for  solving  quasi-one-dimensional  nozzle  flow  are  'resented.  Two  types 
of  characteristic  boundary  conditions,  those  of  i  >er  and  McKenna, 
were  tested  for  each  scheme.  Solutions  of  isentr  :  supersonic  flow 
and  flow  with  shocks  were  obtained  for  a  divergin  zle.  The 
behavior  of  each  boundary  condition  on  both  implicit  and  explicit 
schemes  were  the  same.  They  deviated  from  the  theoretical  values  by 
less  than  one  percent.  A  test  to  determine  the  utility  of  each  scheme 
was  run  by  allowing  the  exit  boundary  conditions  to  change  in  the  hope 
of  forcing  the  shock  to  move  upstream  or  downstream.  The  shock  would 
not  move  in  the  implicit  scheme  for  either  boundary  condition.  The 
explicit  scheme  could  move  the  shock,  but  only  when  Steger's  boundary 
condition  was  used.  (This  is  the  one  which  specified  only  pressure. ) 


I.  INTRODUCTION 


During  the  last  decade  much  progress  has  been  made  both  in  the 
development  of  computer  systems  that  are  more  powerful  and  reliable 
and  in  numerical  techniques  of  solving  the  compressible  Navier-Stokes 
equations.  In  spite  of  this  rapid  progress  in  these  two  related  areas, 
we  still  cannot  calculate  the  flow  fields  past  complete  aircraft 
configurations  (e.g.,  commercial  transport,  fighter  aircraft,  reentry 
vehicles,  etc.)  at  flight  Reynolds  numbers.  If  we  could  efficiently 
calculate  the  flow  field,  there  would  be  no  need  to  conduct  wind 
tunnel  testing  or  build  expensive  experimental  test  vehicles.  These 
costly  experimental  devices  could  all  be  supplanted  by  simply  solving 
the  compressible  Navier-Stokes  equations.  Once  the  technique  of  solv¬ 
ing  these  equations  is  perfected,  the  result  wilx  be  a  dramatic 
decrease  in  the  design  cost  of  new  aerodynamic  vehicles.  Since  we  are 
not  yet  able  to  efficiently  solve  the  compressible  Navier-Stokes 
equations,  continued  research  in  this  area  is  important.  The  process 
of  developing  bigger  and  more  powerful  computers  is  one  aspect  in 
solving  the  Navier-Stokes  equations,  but  the  development  of  numerical 
algorithms  that  efficiently  and  accurately  solve  these  equations  is  the 
motivation  behind  this  paper. 

To  date  a  number  of  algorithms  have  been  written  that  accurately 
solve  for  some  flow  fields.  These  algorithms  however,  are  only  accurate 
for  simplified  configurations  or  components  of  complicated  configura¬ 
tions.  Three  important  factors  that  must  be  considered  in  developing 


a  numerical  algorithm  are:  first,  its  convergence  rate  (how  fast  it 
can  solve  a  problem);  second,  its  robustness  or  reliability  which  means 
providing  correct  answers  for  a  variety  of  problems;  third,  its  ability 
to  adapt  to  complicated  geometries  in  either  two  or  three  dimensions. 

This  criteria  was  and  still  is  the  guide  for  developing  a  new  and  better 
numerical  algorithm. 

Basically,  algorithms  can  be  classified  into  two  categories: 
explicit  and  implicit  methods.  During  the  1950's  and  1960's  several 
explicit  finite-difference  algorithms  were  developed  such  as 
Lax-Wendroff^  types  and  the  popular  MacCormack's  method.^  These 
methods  were  very  popular  and  capable  of  handling  a  variety  of  problems. 
They  were  able  to  solve  high  Reynolds  number  flow  problems  as  well  as 
inviscid-viscous  interactions.  During  the  1970's  however,  researchers 
developed  implicit  methods  to  gain  better  computational  efficiency  over 
the  explicit  methods.  Again,  both  the  implicit  and  explicit  methods 
have  enjoyed  some  degree  of  success,  but  each  have  their  own  unique 
limitations.  For  example,  the  convergence  rate  of  explicit  methods  are 
highly  sensitive  to  space  mesh  size  and  in  order  to  improve  the  solution, 
the  space  meshes  must  be  refined  in  areas  such  as  boundary  surfaces, 
shocks,  and  stagnation  points.  Consequently,  the  computer  time  needed 
to  converge  to  a  solution  becomes  very  large.  Also,  explicit  methods 
are  conditionally  stable.  Therefore,  in  order  to  calculate  high  Reynolds 
number  problems  a  large  amount  of  computing  time  is  required  due  to  the 
small  time  step  limitation  imposed  by  the  stability  condition.  On  the 
other  hand,  implicit  methods  are  unconditionally  stable  and  are  not 
restricted  to  explicit  stability  conditions.  However,  they  are  more 
complicated  than  explicit  methods  and  are  not  as  robust.  They  also 
have  problems  in  handling  the  inviscid  region  as  well  as  with  shock 
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capturing.  Today  the  major  emphasis  is  with  implicit  methods  because 

they  do  not  requive  as  much  computer  time  as  explicit  methods. 

2  3 

(MacCormack  and  Lomax  as  well  as  Hollanders  and  Viviand  have 

presented  papers  that  compare  the  relative  pros  and  cons  of  the  explicit 

and  implicit  methods  as  well  as  their  many  variations.) 

The  purpose  of  this  paper  is  to  study  implicit  techniques  used  in 

solving  the  compressible  Euler  equations.  Areas  such  as  boundary 

conditions,  differencing  techniques  and  shock  capture  are  addressed. 

Using  the  criteria  of  convergence  rate,  robustness,  and  adaptability, 

the  advantages  and  disadvantages  of  implicit  methods  will  be  determined. 


II .  GOVERNING  EQUATIONS 


Throughout  this  paper  the  test  equations  used  are  the  quasi-one- 
ditnensional  Euler  equations.  The  one  dimensional  Euler  equations  are 
used  because  of  their  simplicity  to  demonstrate  the  principles  presented. 
The  extension  to  two  and  three  dimensions  as  well  as  the  full 
compressible  Navier-Stokes  equations  is  straightforward.  It  is 
important  to  stress  the  fact  that  if  the  extension  to  three  dimension 
were  not  straightforward,  then  that  particular  method  would  be  useless. 
Useless,  because  it  cannot  be  applied  to  realistic  problems. 

The  quasi-one-dimensional  unsteady  compressible  form  of  the  Navier- 
Stokes  equations,  neglecting  body  forces  and  heat  sources,  can  be 
written  in  conservation  form  as 


where 


and 


H  = 


0 


(1) 


(2) 


(3) 


and  where 


a 


.3u  ,  3u 

p  "  al 


but  by  using  Stoke's  hypothesis  (X  < 

0  ■  P 


4  3u 
3U  3x 
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(6) 


The  vector  F  is  the  flux  vector  while  U  is  the  vector  of  conservative 
variables  and  H  is  the  source  term.  The  governing  equations  are  in  the 
"weak"  conservation  form  because  the  source  term  (H)  is  outside  of  the 
derivatives  of  the  conserved  variables.  The  primitive  variables  are 
density  p,  pressure  P,  and  velocity  u.  The  total  energy  per  unit  volume 
is  represented  by  e  and  defined  as 


e 


(7) 


where  gamma  (y)  is  the  ratio  of  specific  heats.  The  viscosity 
coefficients  are  X  and  y  while  the  coefficient  of  heat  conductivity 
and  temperature  are  represented  as  k  and  T  respectively. 


Normalized  Governing  Equations 

The  governing  equations  are  normalized  so  that  the  characteristic 
parameters  of  the  problem  may  be  independently  varied.  The  reference 


condition  used  is  the  freestream  condition.  The  dimensional  values 


are  normalized  as  follows 


u 

u  -  — 

u„ 


x 

*  ■  L 


t-rr- 

L/u0 


K  PoUo 


The  continuity  equation  will  be  normalized  to  demonstrate  how  the 
procedure  is  performed. 

Given 


multiply  the  whole  expression  by  — —  to  get 

I. 


3  P  +  3  p  u 
3t  Po  3*  Pou0 
L/uq  L 


or  in  nondimens ional  form 


The  momentum  and  energy  equation  are  nondimensionalized  the  same  way 

only  they  are  multiplied  by  ■  -  y  and  ~~r  respectively. 

Pdun  P0uo 

L  L 
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III.  BOUNDARY  CONDITIONS 


The  study  of  boundary  conditions  is  important  in  solving  any 
partial  differential  equation.  Even  the  extension  to  finite  difference 
equations  does  not  deter  from  this  fact.  The  way  boundary  conditions 
are  chosen  can  influence  the  result  of  any  numerical  scheme  (i.e.  con¬ 
vergence,  stability,  and  accuracy).  In  fact,  though  one  might  start 
out  with  a  stable  interior  scheme  (i.e.,  scheme  for  the  interior  points), 
it  can  become  unstable  and  inaccurate  if  the  boundary  conditions  are 
improperly  treated. 

In  this  section  numerical  boundary  schemes  are  studied. 

Specifically,  those  schemes  that  can  be  used  with  an  implicit  finite 
difference  scheme.  Implicit  boundary  condition  procedures  have  been 
studied  for  several  years  for  the  purpose  of  reaching  a  time-asymptotic 
steady  state  faster  as  well  as  for  allowing  a  time  step  that  is 
unrestricted  by  the  Courant,  Friedrich,  Lewy  (CFL)  stability  criteria. 

Background  Information 

Consider  the  simple  wave  equation 

|^  +  e  |-  -  0,  0Sx<l,  ttf>  (11) 

with  the  initial  condition 

u(x,0)  -  f(x) 

A 

Kreiss  has  shown  that  according  to  mathematical  laws,  in  order  to  be 
able  to  solve  this  type  of  problem,  certain  analytic  boundary  conditions 


must  be  specified.  These  boundary  conditions  are  dependent  upon  the 
sign  of  C  (positive  or  negative).  For  example,  if  C  is  positive  then 
there  must  be  a  boundary  specified  at  X  ■  0,  or  if  C  is  negative  then 
the  boundary  condition  must  be  specified  at  X  *=  1.  In  mathematical 
terms 


u(l,t)  =  g^(t)  for  c<0 
or 

u(0,t)  o  g2(t)  for  c>0 

Once  this  information  is  specified  the  solution  of  this  problem  can 
proceed  in  a  straightforward  manner. 

Now  look  at  a  system  of  two  linear  wave  equations: 


and 
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(12) 


The  two  boundary  conditions  needed  to  solve  this  problem  are  again 
determined  by  the  signs  of  *c'  and  'a'.  For  c>0  and  a>0  the  boundary 
conditions  are  given  by  specifying  u  and  v  at  x  -  0.  Problems  arise 
though  when  'a'  and  'c'  have  different  signs.  To  demonstrate  this 
suppose  c>0  and  a<0.  This  means  that  the  waves  in  u  travel  from  left  to 
right  while  the  waves  in  v  travel  from  right  to  left.  Therefore,  u 
must  be  specified  at  x  -  0  for  one  boundary  condition  while  v  must  be 
specified  at  x  =  1  for  the  second  boundary  condition.  Now  if  u  and  v 
are  expressed  (at  the  boundaries)  in  terms  of  each  other,  problems 
arise.  For  example,  let 


u(x,0)  -  1 


v(x,0)  »  0 


(initial  condition) 


ft 


i 


As  the  right  running  wave  *u'  approaches  the  right  boundary  (x  *  1)  it 
causes  'v*  (which  is  in  terms  of  u)  to  be  non-zero.  A  disturbance  is 
generated  each  time  a  wave  reaches  the  boundary.  This  type  of  boundary 
condition  is  referred  to  as  a  "reflective  boundary  condition". 

The  examples  presented  are  simple  models.  It  was  easy  to  identify 
the  variables  and  where  they  were  to  be  specified.  The  extension  to 
three  coupled  equations  (i.e.  Navier-Stokes)  increases  the  complexity. 
The  choice  of  variables  and  how  they  are  specified  is  no  longer  clear. 
Problems  with  reflective  boundary  conditions  are  another  problem  area. 

A  method  is  needed  to  determine  how  these  complicated  boundary 
conditions  should  be  handled.  This  method  should  represent  the  physics 
and  be  mathematically  correct. 


Characteristic  Theory 

The  mathematical  theory  of  characteristics  for  hyperbolic  systems 
of  equations  is  an  important  clue  in  determining  how  boundary  conditions 
should  be  constructed.  Many  researchers  have  also  used  this  clue 
(ref.  5-9).  Characteristic  theory  is  an  important  clue  because  it  can 
be  used  to  determine  the  number  of  boundary  conditions  needed  to  solve 
a  problem  without  over-determining  it.  Lets  look  again  at  a  model 
hyperbolic  equation 


3u  .  3u 

3t  +  c  51 


o 


(13) 


and  let 


and 
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From  characteristic  theory  it  is  known  that  the  eigenvalues  of  a 
hyperbolic  system  of  equations  (known  as  the  characteristic  direction) 
are  used  to  determine  which  direction  the  information  propagates  at  the 
boundaries.  Knowing  which  direction  the  information  flows  is  vital  in 
knowing  how  many  or  which  boundary  conditions  are  needed.  For  the  model 
equation  the  eigenvalues  (X)  of  the  'C'  matrix  are  determined  by 


or 


X2-l 


0 


finally 


X  ■  ±1 


The  sign  (positive  or  negative)  of  the  eigenvalues  is  what  determines 
which  direction  the  information  flows  (see  Figure  1). 


From  this  example  it  is  clearly  demonstrated  (as  stated  previously) 
that  for  positive  X  which  is  associated  with  r,  a  boundary  condition  is 
needed  at  x  =  0  and  for  negative  X  associated  with  s,  a  boundary 
condition  is  needed  at  x  -  1.  Another  view  of  characteristic  theory  is 
obtained  by  exploring  compatibility  conditions.  More  specifically,  a 
hyperbolic  system  of  equations  such  as  the  one-dimensional  compressible 
Navier-Stokes  equations  is  equivalent  to  three  characteristic 
compatibility  conditions.  These  compatibility  conditions  are  valid 
along  the  characteristic  curves  and  are  in  the  form  of  ordinary 
differential  equations.  The  slopes  of  the  characteristic  curves  are 
given  by  the  eigenvalues.  The  direction  of  these  characteristic  curves, 
which  again  are  determined  by  the  sign  of  the  eigenvalues,  are  used  in 
determining  boundary  conditions.  For  example,  those  curves  that  reach 
the  boundaries  from  inside  the  computational  domain  are  considered 
"admissible".  Admissible  because  they  allow  information  to  propagate 
out  of  the  domain.  Those  curves  that  reach  the  boundary  from  outside 
the  computational  domain  are  called  "inadmissible"  because  information 
cannot  propagate  into  the  domain.  Computational  boundary  conditions 
are  boundary  conditions  that  are  allowed  to  float  or  are  calculated. 
They  are  used  with  compatibility  conditions  associated  with  "admissible" 
characteristic  curves.  Specified  boundary  conditions  are  boundary 
conditions  that  are  specified.  They  are  used  with  those  compatibility 
conditions  associated  with  "inadmissible"  curves.  In  order  to  under¬ 
stand  this  concept  more  easily,  look  at  Figure  2  and  Figure  3. 


outflow  boundary 

I 


outflow  boundary 


points  points  points  points 


Fig.  2  -  Flow  of  Information  Fig.  3  -  Flow  of  Information 

In  Figure  2  the  characteristic  curves  of  A,  B,  and  C  are  all 
admissible,  they  propagate  information  out  of  the  computational  domain. 
Therefore  the  boundary  conditions  must  be  calculated  for  all  three 
compatibility  conditions.  The  sign  of  the  eigenvalues  are  all  positive. 
In  Figure  3  only  curve  C  is  inadmissible,  therefore  a  specified 
boundary  condition  is  needed  for  this  compatibility  condition  while  the 
boundary  conditions  for  curves  A  and  B  are  again  calculated.  Here,  A 
is  positive  for  curves  A  and  B  but  negative  for  the  'C"  curve.  It  is 
important  to  note  that  the  eigenvalues  do  not  need  to  be  calculated. 

The  only  information  needed  is  the  sign  (positive  or  negative).  The 
role  of  positive  and  negative  eigenvalues  at  the  left  boundary  are 
reversed  at  the  right  boundary.  These  concepts  can  and  will  be  used 
with  the  Navier-Stokes  equations. 

Nozzle  Boundary  Conditions 

The  determination  of  boundary  conditions  for  a  nozzle  is  dictated 
by  characteristic  theory.  The  slope  of  the  characteristic  curve  or 
sign  of  the  eigenvalue  determines  what  variables  are  specified  and 


which  are  calculated.  As  will  be  shown  later,  the  eigenvalues  for  the 
one-dimensional  inviscid  compressible  Navier-Stokes  equation  are: 

A^  ®  u 
A2  -  u+c 
A^  ■  u-c 

where  "u"  is  the  velocity  and  "c"  is  the  speed  of  sound.  It  is  evident 
that  the  eigenvalues,  which  determine  the  boundary  conditions,  will  be 
determined  by  the  velocity  (i.e.,  subsonic  or  supersonic). 

For  nozzle  flow,  the  boundary  conditions  are  simply  determined  by 
the  velocities  at  the  inflow  and  at  the  exit.  For  the  subsonic  inflow 
condition  the  signs  of  the  eigenvalues  are: 

A^  ■  u>0 
A^  ■  u+c>0 
Aj  «  u-c<0 

here  A^  and  A2  are  positive  while  A^  is  negative.  This  corresponds 
to  the  slopes  of  Figure  3  (only  this  time  it  is  for  the  inflow  boundary) 
The  only  characteristic  curve  that  propagates  information  from  the 
computational  domain  to  the  boundary  is  \y  This  means  that  the 
variables  associated  with  A^  should  be  calculated  while  the  variables 
associated  with  A^  and  A2  should  be  specified.  The  variable  associated 
with  A^  is  from  the  continuity  equation.  The  variables  associated  with 
A2  and  A^  are  derived  from  the  momentum  and  energy  equations 
respectively.  For  the  supersonic  inflow  condition: 

^  ■  u>0 
^5  ■  u+c>0 


all  the  eigenvalues  are  positive  which  correspond  to  characteristic 
curves  propagating  information  to  the  boundary  from  outside  the 
computational  domain.  Therefore  all  the  inflow  variables  should  be 
specified.  Note  that  the  eigenvalues  for  the  subsonic  outflow  have  the 
same  sign  as  the  eigenvalues  for  the  subsonic  inflow  but  the  boundary 
conditions  are  not  the  same.  They  are  reversed.  For  the  subsonic 
outflow  according  to  Figure  3  X^  and  should  be  calculated  while  X^ 
should  be  specified.  For  supersonic  outflow  conditions  all  the 
variables  should  be  calculated.  This  was  the  procedure  used  in 
determining  boundary  conditions  for  the  algorithm. 

Note,  nothing  has  been  said  about  the  variables  associated  with  the 
eigenvalues.  At  first,  one  would  think  that  the  variables  should  be 
the  primitive  variables:  density-continuity  equation,  velocity-momentum 
equation,  and  pressure-energy  equation.  Upon  closer  examination  it  is 
noticed  that  when  the  equations  (Navier-Stokes)  are  uncoupled  these  are 
not  the  correct  variables. 


Characteristic  Variables 

Characteristic  variables  are  the  variables  associated  with 
eigenvalues.  These  are  the  variables  that  need  to  be  specified  or 
calculated.  The  inviscid  Navier-Stokes  equation  representative  of  one¬ 
dimensional  flow  is  given  as 


3U  3F 

it  *■& 


+ 


H  = 
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(14) 


Now  by  using  the  Jacobian  Matrix  which  is  defined  as 


$((Y-1)  ~  uc)  B(c  -  (y~1)u) 
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0((Y-D  ~2  +  uc)  -0(c  +  (Y“l)u) 
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Now  substitute  equation  17  into  equation  16  to  get 

3u  .  3(xDx_1)u  A 


ft  + 


3x 


-1 


Now  multiply  the  whole  equation (22)by  X  to  get 

-1  3u  .  -1  3(xDx  l)n  n 

*  5t  ♦*  - 35 - 0 

Since  X  *  is  a  constant  matrix  (locally)  it  can  be  placed  inside  the 
derivative  expression. 


or 


where 


and 


3x  *u  3Dx  *u 
3t  3x 


X"1  U  -  Q 
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This  equation 


looks  like  the  simple  wave  equation  (12) that  was  used  as  a  model. 
Instead  of  the  variables  u  and  v  the  variables  are  now  Q^,  Qj,  and  Q^, 
which  is  simply  the  U  matrix  (vector  of  conservative  variables) 
multiplied  by  the  eigenvector  matrix  X  The  results  are 


(29) 


It  was  assumed  that  deviations  from  the  free  stream  are  so  small  that 


the  entries  in  the  Jacobian  can  be  treated  as  constants-locally .  That 
is  where  the  O-subscript  originates.  The  eigenvalues  associated  with 
'a'  and  'c'  are  u,  u+c,  and  u-c. 

From  the  sign  of  the  eigenvalues  (u,  u+c,  u-c)  and  the 
characteristic  variables,  boundary  conditions  can  be  defined. 


Numerical  Conditions 

Numerical  boundary  conditions  are  boundary  conditions  used  with 
finite  difference  equations.  They  are  the  numerical  implementation  of 
analytic  boundary  conditions.  Interesting  problems  arise  between  the 
differencing  technique  and  the  way  boundary  conditions  are  chosen.  For 
example,  the  way  the  spatial  derivative  term  is  differenced 

affects  the  way  boundary  conditions  are  specified.  To  demonstrate 
this  fact,  suppose  a  central  finite  difference  approximation  is  used 
for  the  spatial  term.  There  are  now  two  points  outside  the  boundaries 


that  need  to  be  either  specified  or  calculated.  Another  way  to  look 
at  it  is  to  look  at  a  discrete  mesh  (Figure  4) 


inflow  boundary  outflow  boundary 


-1  0  1  2  3  4  . J-l  J  J+l 


Fig.  4  -  Discrete  Computational  Mesh 


with  J  points  in  the  computational  domain.  The  points  -1  (outside  the 

inflow  boundary)  and  J+l  (outside  the  outflow  boundary)  are  outside  the 

computational  domain.  The  value  of  these  points  must  be  made  known 

til  t  J) 

in  order  to  solve  for  the  J  and  0L  point  according  to  the  difference 
relation 


3u  ("j.l  -  “j-l> 

15x  j  2  Ax 


(30) 


and 


(31) 


How  these  values  are  determined  is  directly  related  to  the  boundary 
conditions. 

As  already  stated  before,  the  sign  of  the  eigenvalues  determines 
which  direction  information  flows  at  the  boundary.  Boundary  points  are 
calculated  by  use  of  this  information.  For  example,  if  the  eigenvalues 
are  positive  at  the  outflow  boundary,  a  backward  differencing  scheme 
should  be  used.  If  the  eigenvalues  are  negative  then  a  forward 
differencing  scheme  should  be  used.  This  is  because  with  positive 
eigenvalues  information  flows  from  inside  the  computational  domain  to 
the  boundary.  Therefore  one  would  not  want  to  use  a  forward 
differencing  scheme  because  then  there  would  not  be  any  continuity 


interior  points  exterior  points  interior  points  exterior  points 


X>0 


A<0 


backward  differencing  scheme 


forward  differencing  scheme 


Fig.  5  -  Differencing  Scheme 


IV.  MACCORMACK'S  IMPLICIT  ALGORITHM 


MacCormack  has  developed  a  new  implicit  algorithm  for  solving  the  inviscid 
Navier-Stokes  equation.  It  is  an  extension  of  his  1969  explicit 
predictor-corrector  algorithm^.  Because  it  is  an  extension  of  an 
explicit  algorithm  this  new  scheme  is  not  truly  implicit.  For  example, 
it  is  not  unconditionally  stable  but  it  does  retain  a  few  of  the 
advantages  of  implicit  algorithms  such  as  larger  time  steps. 

The  new  method^  is  composed  of  two  parts.  The  first  part 
explicitly  solves  the  governing  equations  by  using  known  flowfield 
properties  to  determine  local  flowfield  changes.  The  second  step 
numerically  transforms  the  equation  into  an  implicit  form  to  remove 
the  explicit  stability  criteria  (At  ■  ^u+'a]  )*  7116  true  flowfield 

changes  are  then  calculated  and  substituted  into  a  Taylor  series  in 
time  to  find  the  new  flowfield  properties. 


Derivation 


Differentiating  the  one -dimensional  Navier-Stokes  equations  with 
respect  to  time  gives 


Using  the  Jacobian,  the  above  equation  can  be  rewritten  as 


(32) 


(33) 


A  closer  look  at  equation  (35)  reveals  the  two  step  procedure. 
The  right  hand  term  represents  the  explicit  solution  of  the  governing 
equations  while  the  left  hand  side  of  the  equation  represents  the 
numerics  of  the  implicit  algorithm.  Because  MacCormack  uses  a 
predictor-corrector  approach  equation  (35)  is  rewritten  as 
predictor: 
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corrector: 


and  where 


3H 

w 


0  0  0 

(y-1)  ..2  3A  (Y-l)  .  3A  (Y-l)  3A 
2A~  U  17  "  "A  u  17  "A 


0 
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Spatial  derivative  terms  are  approximated  by  a  forward  difference  for 
predictor  and  backward  difference  for  corrector. 

The  Jacobian  (A)  is  the  same  as  previously  defined  in  equation  17 
and  18.  The  absolute  value  of  the  Jacobian  (|a|)  which  is  used  in  the 
predictor  and  corrector  equations  is  defined  as 

|A|  *  X  A a  X-1  (38) 


(The  transformation  matrices  have  already  been  given,  see  equations 
20  and  21). 

The  eigenvalue  matrix  is  written  as 
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0 
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(39) 
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where 


\  •  (l“l  (If)  .  ».o) 

>A2  ■  ■“  (lu'fol  "CE(ff)'  0,°) 

\  ■  Bax  (l u"c i  -CE(H).  o-o) 

CE  =  a  constant  related  to  the  Courant  number  used  for 
the  explicit  stability. 


where 


At  < 


CE  (Ax) 
(u+c) 


If  the  explicit  stability  criteria  is  used  the  X^  matrix  goes  to  zero. 
This  allows  the  true  flowfield  changes  to  be  calculated  in  the  explicit 
part  of  the  program  thereby  skipping  altogether  the  implicit  part. 

Throughout  this  paper  only  the  inviscid  form  of  the  equations  were 
used.  If  viscous  effects  are  included,  the  eigenvalue  matrix  would 
become: 
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(Note:  v  is  nondimensionalized  by  dividing  it  by  v0.) 

Modification  to  MacCormack's  Implicit  Algorithm 

Careful  analysis  of  MacCormack's  predictor  -  corrector  equations 
reveals  that  after  the  explicit  (local)  flowfield  changes  have  been 
calculated  and  substituted  into  the  right  hand  side  of  equation  (35) 
an  upper  block  bidiagonal  system  of  equations  will  result  for  the 


23 


predictor  step.  A  lower  block  bidiagonal  system  of  equations  will 
result  for  the  corrector  step.  In  order  for  the  solution  to  progress  a 
matrix  inversion  must  be  performed  at  each  computational  point  which 
will  result  in  a  large  amount  of  computer  time,  thus  defeating  one 
advantage  of  implicit  methods. 

MacCormack  however,  has  devised  an  ingenious  way  to  circumvent 

the  inversion  process.  The  [ A (  matrix,  which  needs  to  be  inverted,  is 

diagonalized  thus  making  its  inversion  trivial.  Problems  with  this 

technique  arise  when  the  source  term  "H"  is  included  in  the  implicitly 

approximated  governing  equations  (equation  35).  The  diagonalization 

cannot  be  performed  and  a  matrix  inversion  must  be  computed  for  each 

point.  However,  a  modification  to  MacCormack's  scheme  that  will  allow 

a  diagonalization  by  neglecting  the  source  term  in  the  implicit 

12 

numerics  was  presented  in  a  paper  by  White  and  Anderson  .  They 
reasoned  that  the  source  term  "H"  could  be  removed  from  the  implicit 
numerics  because  "The  physics  of  the  flow  is  carried  by  the  explicit 
step  (which  contains  the  source  term  in  the  governing  equations),  and 
the  implicit  step  (which  under  the  present  modification  dees  not  contain 
the  source  term)  is  simply  numerics  to  bring  about  enhanced  stability." 
Using  their  modification  the  implicitly  approximated  equation  can  be 
rewritten  as  predictor: 


l  ii 


24 


corrector: 
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This  is  in  a  form  that  results  in  an  upper  and  lower  block 
bidiagonal  system  which  will  allow  a  diagonalization  to  occur. 


Solution  Method 

As  already  stated  before,  this  scheme  is  a  predictor  -  corrector 
type.  The  predictor  equations  are  solved  first  and  their  results 
are  used  in  solving  the  corrector  equations.  The  results  from  the 
corrector  step,  the  true  flowfield  properties,  are  then  used  in 
calculating  the  predictor  step.  This  process  is  repeated  until  the 
scheme  converges.  In  order  to  understand  this  scheme  more  fully  the 
predictor  step  will  be  reviewed  in  detail.  The  principles  used  in  the 
corrector  step  are  the  same  as  those  used  in  the  predictor  step.  The 
extension  therefore  should  be  straightforward. 

Initially  all  the  conservative  variables  "U"  are  specified 
(initial  guess)  for  all  i  «  1,2,3...J  (J  *=  is  the  total  number  of 
computational  points).  Once  the  end  boundary  condition  is  found 
| A |  the  right  hand  side  of  the  predictor  equation  can  be 

solved  and  represented  as  "W" 


W  ■  AU*?  +4”  |AI  ^  i 

1  Ax  1  1  l+l  i+i 


(41) 
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In  order  to  solve  for  6U^+*,  a  matrix  inversion  of  (i  +  ^  ( A [  ^ )  must 
be  performed.  By  using  the  relation  A  =  XDX  ‘  the  matrix  may  be 
diagonalized  thereby  foregoing  the  costly  (in  computer  time)  matrix 
inversion.  Making  this  substitution,  the  equations  can  be  written  as 

(!  ♦  S  X  \  X'1)  <l  '  "  <«> 

Now  pre-multiply  both  sides  by  X  *  to  get 


(x-1  +  Da  X-1)  SU^  -  x"1  W  -  V 


(43) 


Regrouping,  the  left  hand  side  is  now  written  as 


(x  *  S  da)  X_1  “T  -  v 


(44) 


The  matrix  ^1  +  ~  is  diagonal  so  that  its  inversion  is  trivial, 
The  equations  can  now  be  written  as 


X-1  6uN+1  =  ^  +  d^~1  v  =  y  (45) 

n+t  . 

Where  <51K  is  found  by  multiplying  the  equation  by  X.  The  flowfield 

changes  for  the  point  i  have  just  been  calculated.  Now,  in  order  to 

calculate  the  values  for  the  next  point  ( i— 1 )  the  term  { A |  6U^+^  must 

be  calculated.  Calculating  this  term  will  give  all  the  values  for  the 

right  hand  side  of  the  equation.  Once  all  these  values  are  obtained  the 

N+l 

procedure  used  to  solve  for  6lL  ^  (on  the  left  hand  side)  is  exactly 

N  N+T 

the  same  as  stated  before.  An  efficient  way  to  calculate  the  | A |  ^  61K 
term  is  to  multiply  the  already  known  value  of  Y  by  then  by  X.  In 
mathematical  terms 


X  D,  Y 
A 


(46) 
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Using  this  method  the  predictor  step  sweeps  down  from  i  «  J  to  the 


i  ■  1  computational  point.  The  corrector  step  starts  at  i  =  1  and 
marches  to  the  i  ■  J  point.  Local  values  are  used  in  calculating  the 
matricies  X,  D,  and  X  *  at  each  point. 

The  solution  process  can  be  summarized  in  the  following  seven  steps 


1. 

2. 

3. 

4. 

5. 

6. 

7. 


M  -  K  *  §  M  L  <1 


V  °  x”1  w 


D.  is  calculated 
A 


Y  -  (I  + 


(l  *  s  da) 

6U*J+1  =  X  Y 

Z  -  D.  Y 
A 


-1 


J A)  ?  6uJ+1  «XZ 


Because  the  matrix  is  diagonal  in  step  4  the  inversion  is  trivial. 
Step  5  calculates  the  flowfield  changes  while  step  7  calculates  the 
flux  to  be  used  for  the  next  computational  point.  If  Neumann  type 
boundary  conditions  are  used  the  [ A |  2  5U2^+^  term  is  saved  and  used 
as  the  boundary  condition  for  the  corrector  step.  The  1*1  r.i  <1 
term  is  saved  from  the  corrector  step  and  used  as  the  boundary 
condition  for  the  predictor  step.  Because  of  this  the  boundary 
conditions  are  not  truly  implicit  but  lag  by  1/2  time-step. 


Extension  to  Multi-Dimensions 

The  value  of  any  one -dimensional  scheme  is  proportional  to  its 
ability  to  adapt  to  multi-dimensions.  A  one-dimensional  scheme  that 
is  extremely  accurate  but  not  able  to  convert  to  multi-dimensions  is 


severely  limited.  The  method  in  devising  any  numerical  scheme  is 
to  start  out  with  a  simple  model,  then  add  upon  that  model  so  that  it 
can  handle  a  wider  variety  of  problems.  If  the  model  cannot  expand  then 
it  is  a  poor  model.  If  the  model  can  expand  while  still  retaining  some 
degree  of  accuracy  then  one  has  a  good  model.  The  idea  is  to  start  out 
simple  and  then  become  more  sophisticated. 

MacCormack's  implicit  scheme  is  a  robust  scheme.  It  can  handle  a 
wide  variety  of  problems  as  well  as  being  able  to  expand  to  two  and 
three  dimensions.  The  extension  from  the  simple  one-dimensional  case 
to  two  and  three  dimensions  is  straightforward.  To  demonstrate  this 
fact  the  scheme  will  be  developed  for  the  two-dimensional  case.  The 
general  equation  form  is  given  by 


3U  A  3F  3G 
It  lx  3y 

Differentiating  the  general  equation  (48) 
using  the  Jacobian  yields 


-  0 

with  respect  to  time  and 


(48) 


(49) 


(A  and  B  are  the  Jacobians).  Now  implicitly  approximate  the  above 
equation  (49)  in  time  to  get 


(‘ 


3  A  • 
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(50) 


The  dots  have  the  same  meanii  as  before  (the  derivative  acts  on  the 
term  outside  the  parenthesis).  Uj.lng  the  following  substitutions 


aun 


At 


au” 

TT 


(51) 
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The  solution  procedure  for  this  scheme  in  two-dimensions  is  explained 
in  detail  in  reference  11.  It  will  not  be  discussed  further. 


From  this  example  it  has  been  shown  that  MacCormack's  implicit 
scheme  can  be  expanded  to  multi-dimensions.  The  extension  to  three- 
dimensions  is  straightforward.  This  little  example  shows  the  utility 
of  this  scheme. 

Boundary  Conditions 

The  boundary  conditions  used  in  MacCormack's  algorithm  are  based 

on  characteristic  theory.  They  are  applied  to  a  diverging  nozzle  that 

has  a  constant  supersonic  inflow  and  an  outflow  that  varies  between 

13 

subsonic  and  supersonic.  Two  boundary  conditions,  one  by  McKenna  and 

14 

the  other  by  Steger  ,  were  studied  for  their  use  in  subsonic  (exit) 
flows . 

The  upstream  boundary  conditions  for  the  supersonic  (inflow) 
nozzle  were  set  by  holding  the  characteristic  variables  constant.  This 
is  because  the  eigenvalues  which  are  all  positive  dictate  specified 
boundary  conditions.  The  characteristic  variables  which  are  functions 
of  the  primitive  variables  (density,  velocity,  and  pressure)  are  the 


variables  that  need  to  be  specified.  They  are  specified  by  holding  the 
primitive  variables  constant  at  their  inflow  values.  They  are 
represented  as 


inflow 

mathematically  numerically 

P 

<?1  =  p  -  •jjr-r  =  constant  p(l)  =  constant 

=  P  +  P0C0U  =  constant  U(l)  =  constant 

Q3  =  P  -  p0C0U  =  constant  PCD  =  constant 

Numerically  these  boundary  conditions  were  written  so  that  the  inflow 
variables  were  the  same  for  each  time-step. 

The  downstream  boundary  conditions  for  the  supersonic  (exit) 
nozzle  were  calculated.  The  eigenvalues  are  positive,  thereby 
indicating  right  running  waves,  the  same  as  at  the  inlet.  The  difference 
is  that  the  right  running  waves  approach  the  boundary  from  inside  the 
computational  domain.  The  characteristic  variables  must  therefore  be 
calculated.  Neumann  type  boundary  conditions  were  used  for  this 
condition.  They  are  represented  as 

outflow 


mathematically 

numerically 
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Using  this  type  of  boundary  condition  insures  that  waves  are  not 
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reflected  at  the  boundary  . 

The  exit  boundary  conditions  (Neumann)  were  implemented  into  the 
program  by  setting  the  next  to  last  computational  point  in  the  corrector 
step  equal  to  the  last  point  in  the  predictor  step. 

(Vw  ■  <Vw  - 1 

Predictor  Corrector 


This  is  because  the  predictor  step  sweeps  from  the  exit  of  the  nozzle  to 
the  entrance  while  the  corrector  sweeps  from  the  entrance  to  the  exit. 
Therefore,  the  exit  boundary  condition  (used  in  the  predictor  step) 
should  come  from  the  corrector  step. 

The  boundary  conditions  of  McKenna  and  Steger  were  used  for  the 
subsonic  (exit)  boundary  condition.  The  subsonic  boundary  condition  is 
more  complicated  because  two  of  the  eigenvalues  are  positive  while  the 
third  is  negative.  This  corresponds  to  having  the  first  two  variables 
calculated  while  the  third  is  specified.  Mathematically  it  is 
represented  as 


3Q1 

IT  “ 


o 


0 


n  »  Constant  =  K3 


Numerically,  in  order  to  specify  K3  which  is  equal  to 


K3  -  P  -  p0C0U 


the  pressure  and  velocity  must  be  known  at  the  outflow  boundary. 
(Remember  p0  and  C0  are  values  used  from  the  previous  time  step.) 

This  is  not  helpful  for  real  nozzle  flow  because  the  only  known 
variable  downstream  is  the  pressure.  Specifying  the  velocity  downstream 
is  mathematically  correct  but  not  realistic.  McKenna  has  used  this 
information  in  writing  his  boundary  conditions.  They  are  given  as 


n  Do  /  3  j-l  .  „n 

pj-l  +  2Cq"  I  p^ST  -  PoCT  +  Vl 


0^  =  7  |U;  ,  + 


j  2  lj-1  p0c0  A.p0c0y 

'  /pn  K 

P?  „  +  u1}  ,  +  _L_ 

j  2  l  p0C0  j-l  AjP0C0 


where 


K3  -  P.  -  P0C0U. 

J  J 

(The  derivation  of  this  boundary  condition  is  in  the  Appendix) 

One  assumption  that  McKenna  makes  that  is  suspect  is  his  assumption 
that  K3  is  constant. 

A  set  of  subsonic  downstream  boundary  conditions  which  only 
specify  the  pressure  were  also  used.  They  are  Steger’s  boundary 
conditions.  The  derivations  of  these  boundary  conditions  are  also 
given  in  the  appendix.  They  are  represented  as 


3.  *  p.  .  +  7T  (Pco  ■  P-  ,) 


vi  ’  Vi  *  577 ( Vi  -  p~> 


p,  -  p« 


This  type  of  boundary  condition  is  more  common  in  an  experimental  setup 


but  allows  reflection  of  waves 


V.  RESULTS 


A  series  of  computer  runs  were  made  using  both  MacCormack’s 
implicit  and  explicit  algorithms.  Boundary  conditions  were  tested  for 
various  nozzle  conditions  to  determine  their  impact  on  both  schemes.  A 
problem  of  special  interest  that  was  investigated  was  to  determine  how 
MacCormack's  algorithm  (with  the  help  of  characteristic  boundary 
conditions)  could  handle  the  relationship  between  the  shock  location 
and  back  pressure. 

The  results  are  plotted  against  the  theoretical  solution  for 
comparison.  Notice  that  the  theoretical  shock  is  smeared.  This  is 
because  the  theoretical  solution  is  calculated  and  plotted  at  each 
computational  point.  This  will  give  a  better  comparison  between  the 
numerical  results  and  the  theoretical  solution. 

The  first  test  performed  was  comparing  the  implicit  to  the  explicit 
algorithm  for  the  supersonic  nozzle.  Boundary  conditions  were  not  a 
factor  in  this  computer  run  because  both  algorithms  handled  the 
conditions  the  same — numerically.  The  supersonic  downstream  boundary 
conditions  were  Neumann  type.  The  results  were  not  surprising.  The 
implicit  run  converged  in  135  time-steps  while  the  explicit  converged 
in  1408  time-steps.  Interestingly,  the  implicit  run  did  not  converge 
to  the  same  values  of  the  explicit  run.  This  was  a  result  of  a  constant 
which  is  related  to  the  CFL  condition  located  in  the  eigenvalue  matrix. 
After  this  constant  was  adjusted  to  0.9  both  algorithms  converged  to 
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the  same  values.  These  values  were  within  1%  of  the  theoretical  values. 
The  results  are  plotted  on  Figures  6  through  9. 

The  second  test  performed  was  a  comparison  of  McKenna's  and 
Steger's  boundary  conditions.  McKenna's  was  the  characteristic  boundary 
condition  where  both  pressure  and  velocity  are  specified  while  Steger's 
characteristic  boundary  condition  only  specified  the  pressure.  Each 
boundary  condition  was  used  in  the  implicit  and  explicit  algorithm. 

The  subsonic  (exit)  nozzle  was  used  because  characteristic  theory  states 
that  with  subsonic  flow  there  will  be  one  left  running  wave  and  two 
right  running  waves.  This  is  a  perfect  condition  in  which  to  test  the 
two  boundary  conditions  for  their  reflection  effects  on  the  algorithm. 
The  convergence  of  both  boundary  conditions  in  the  implicit  case  compare 
quite  well.  The  boundary  condition  that  specified  pressure  and  velocity 
(McKenna)  required  one  less  time-step  to  converge  than  Steger's  boundary 
condition.  Except  for  the  jumps  (Gibb's  phenomena)  around  the  shock, 
both  conditions  converged  to  the  theoretical  values  with  a  deviation  of 
less  than  one  percent.  The  results  are  plotted  in  Figures  10  through  13 
The  plots  for  the  explicit  case  are  located  in  Figures  14  through  17. 
Notice  that  the  Figures  are  not  as  smooth  as  those  of  the  implicit  case. 
After  2,500  iterations  the  scheme  still  had  not  converged.  It  simply 
oscillated  around  the  theoretical  values.  It  was  suggested  that  by 
allowing  the  time-step  At  to  float  it  would  cause  the  waves  to 
decrease.  Accordingly,  At  was  set  to 


At 


CN  »  AX 
U  +  C 


for  each  computational  point.  The  results  justified  this  line  of 
reasoning.  Note  how  smooth  the  results  are — just  like  in  the  implicit 
case  (see  Figures  18  and  19).  The  reason  the  jumps  disappeared  is 
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Figure  6.  Mach  Number  Distribution  for  a  Supersonic  Nozzle 
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Figure  7.  pressure  Distribution  for  a  Supersonic  Nozzle 
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Figure  8.  Each  Number  Distribution  for  a  Supersonic  Nozzle 
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Figure  11.  Pressure  Distribution  for  Subsonic  Flow  with  a  Shock 
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rigure  12.  bach  number  Distribution  for  Subsonic  Flow  with  a  shock 
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Figure  13*  Pressure  Distribution  for  Subsonic  Flow  with 
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Figure  15.  Pressure  Distribution  for  Subsonic  Flow  with  a  Shock 
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Figure  16.  i.ach  Number  Distribution  for  Subsonic  Flow  with  a  Shock 
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Figure  1?.  Pressure  Distribution  for  Subsonic  Flow  with  a  Shock 
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Figure  18.  Mach  Number  Distribution  for  Subsonic  Flow  with  a  Shock 
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Figure  19.  Pressure  Distribution  for  Subsonic  Flow  with  a  Shock 


because  the  truncation  error  from  the  time  derivative  term  cancelled  the 
artificial  viscosity  generated  by  the  convective  term.  (For  further 
information  concerning  this  point  see  ref.  15). 

The  results  up  to  this  point  showed  that  the  affects  of  both 
boundary  conditions  are  about  the  same. 

A  special  test  was  run  to  determine  how  MacCormack's  algorithms 
could  handle  the  relationship  between  shock  location  and  back  pressure. 
Both  boundary  conditions  were  used  to  test  their  effects  on  the 
algorithms.  The  results  are  located  in  Figures  20  through  27.  The 
test  was  run  by  allowing  the  scheme  to  converge  to  an  initial  guess. 
Then,  after  convergence  was  achieved  the  boundary  condition  was  changed. 
It  was  initially  hoped  that  the  implicit  scheme  would  allow  the  shock 
to  move  upstream  especially  when  McKenna's  boundary  condition  was  used. 
This  condition  was  the  one  that  specified  pressure  and  velocity  at  the 
exit.  It  was  thought  that  by  specifying  two  variables  at  the  boundary 
it  would  "force"  the  shock  to  move,  but  such  was  not  the  case. 

The  first  test  case  involved  the  implicit  algorithm.  McKenna's 
boundary  condition  did  not  work  for  either  moving  the  shock  upstream  or 
downstream.  Immediately  after  the  initial  guess  converged  and  the 
boundary  condition  changed,  negative  values  for  the  pressure  and  density 
appeared.  The  results  were  meaningless.  The  results  for  the  implicit 
algorithm  with  Steger's  boundary  conditions  were  the  same  as  McKenna's. 
They  both  failed.  They  did  the  same  thing — produced  negative  values 
for  pressure  and  density  and  became  numerically  unstable.  This  was 
surprising.  It  was  hoped  that  MacCormack's  implicit  algorithm  could 
handle  the  problem  of  allowing  the  shock  to  move. 


The  results  for  the  explicit  case  were  successful.  McKenna's 
boundary  condition  that  specified  pressure  and  velocity  still  did  not 
allow  the  shock  to  move  upstream  or  downstream,  but  Steger's  boundary 
condition  that  specified  just  the  pressure  did.  As  Figures  20  through 
23  show,  Steger's  boundary  condition  moves  the  shock  up  and  down  the 
nozzle  right  where  it  is  supposed  to  be  (for  the  corresponding  pressure) 
This  is  an  interesting  point.  It  shows  that  the  boundary  condition 
that  specifies  only  pressure  is  the  logical  choice  to  use  both 
theoretically  and  numerically.  Theoretically,  pressure  is  the  driving 
mechanism  for  nozzle  flow  problems.  Numerically,  specifying  pressure 
at  the  boundary  works . 

It  was  then  decided  to  allow  the  time-step  (At)  to  float  to  see  if 
the  jumps  could  be  smoothed  out.  The  test  was  run  for  both  boundary 
conditions  in  the  explicit  algorithm.  The  results  for  attempting  to 
move  the  shock  upstream  are  located  in  Figures  24  and  25.  Only  Steger's 
boundary  condition  was  stable.  Notice,  the  shock  did  not  move.  In 
fact,  the  results  are  exactly  the  same  from  the  inflow  to  the  shock 
location  (original  location)  as  the  initial  theoretical  values.  Once 
past  the  shock  the  results  for  each  boundary  condition  varied. 

Analyzing  the  results  it  appears  as  though  once  the  shock  initially 
forms,  a  change  in  the  downstream  boundary  condition  cannot  be  felt 
upstream  of  the  shock.  The  effects  of  the  downstream  boundary 
conditions  are  just  reflected  back  and  forth  from  the  shock  to  the  exit. 

The  attempts  to  let  the  shock  move  downstream  were  disappointing. 
Neither  boundary  condition  worked.  The  condition  that  specified  only 
pressure  did  converge  but  again  the  shock  did  not  move.  McKenna's 
boundary  condition  was  again  unstable.  The  upstream  values  were  the 
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Figure  20.  Mach  Number  Distribution  for  Moving  the  Shock  Upstream 
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Figure  21.  Pressure  Distribution  for  Moving  the  Shock  Upstream 
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Figure  22.  Mach 
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Figure  23.  Pressure  Distribution  for  hoving  the  Shock  Downstream 
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Figure  24.  Mach  Number  Distribution  for  Trying  to  Move  the 
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Figure  25.  Pressure  Distribution  for  Trying  to  Move  the  Shock  Upstream 


same  as  the  initial  guess  (converged  solution),  but  the  downstream 
values  were  radically  different  (same  as  before).  See  Figures  26  and  27 
for  the  results.  The  results  for  this  section  were  very  surprising. 

It  was  thought  that  the  shock  would  move  downstream  (maybe  not  upstream, 
but  at  least  downstream)  due  to  the  subsonic  nature  of  the  flow,  but 
as  stated,  the  results  readily  proved  this  inaccurate.  This  is  an 
interesting  case  because  when  At  is  held  constant  the  shock  moves! 
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Figure  26.  Kach  Number  Distribution  for  Trying  to  Move  the  Shock  Downstream 
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VI.  CONCLUSIONS 


Numerical  experiments  for  the  quasi-one-dimensional  Euler 
equations  have  been  applied  to  nozzle  flow  problems.  These  experiments 
included  boundary  condition  analysis  and  a  comparison  of  MacCormack's 
implicit  and  explicit  algorithms.  The  results  from  these  experiments 
are  judged  by  how  well  the  equations  are  solved.  The  judgment  criteria 
is  convergence  rate,  robustness,  and  adaptability. 

There  were  two  (exit)  boundary  conditions  used,  each  based  on 
characteristic  theory.  One  boundary  condition  specified  pressure  and 
velocity  (characteristic  variable)  while  the  other  only  specified 
pressure.  The  convergence  rate  for  both  conditions  was  just  about  the 
same.  The  condition  that  specified  pressure  and  velocity  was  initially 
thought  to  be  mathematically  superior  in  that  a  non-reflecting  boundary 
condition  would  prevail.  The  results  showed  just  the  opposite.  The 
condition  that  only  specified  pressure  was  robust.  This  was  the  only 
boundary  condition  in  which  the  shock  moved  upstream  and  downstream. 

The  adaptability  of  both  boundary  conditions  are  about  the  same.  This 
is  a  result  of  the  fact  that  they  were  both  derived  from  characteristic 
theory.  Judging  the  two  boundary  conditions,  it  is  evident  that  the 
condition  which  specifies  only  pressure  is  more  useful. 

Each  of  MacCormack's  algorithms  has  their  advantages  and 
disadvantages.  The  explicit  method  is  very  reliable  and  can  be  used 
for  a  wide  variety  of  problems.  The  main  drawback  is  the  large  amount 
of  computing  time  required  to  solve  a  problem.  MacCormack's  implicit 


method  is  a  new  method  which  yields  highly  accurate  results  with  a 
considerable  amount  of  reduction  in  comp  iter  time.  The  method  is  not 
truly  implicit  because  it  is  not  unconditionally  stable.  The  main 
advantage  of  the  implicit  method  is  that  when  the  problem  becomes  more 
complex  (viscous  effects  and  turbulence)  the  explicit  part  of  the 
scheme  becomes  more  complex  too,  but  the  implicit  numerics  are  not 
affected — therefore,  a  great  reduction  in  computer  time  exists.  The 
inability  of  the  implicit  scheme  to  move  the  shock  upstream  or 
downstream  was  a  disappointment. 

Areas  of  further  study  that  need  to  be  investigated  are: 

1)  McKenna's  boundary  condition  and  why  they  did  not  work  for  the 
explicit  case.  2)  The  importance  of  initial  conditions.  (This  might 
account  for  the  fact  that  McKenna's  boundary  condition  did  not  work.) 
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APPENDIX  A 


Derivation  of  the  Jacobian 


The  Jacobian  (A)  is  defined  as 


Ai»<u> 


In  order  to  calculate  it,  the  flux  vector  "F"  should  be  redefined  in 
terms  of  the  conservative  variables  "U". 

Given 


U  »  {  pu  i  =  {  m 


then  "F"  can  be  rewritten  as 


F  -  |p  +  pu2  l  -  J  (f-l)e  + 


[u  (e  +  p)J 


Y  me  _  CY-1)  m3 
P  2 


The  Jacobian  can  now  be  computed  as 
A11  3U1  ^  0 


a  1  3®  1 

A12  3U2  "  ‘§5  *  1 


a  1  3m  n 

A13  "  1u7  "  37  “  ° 


APPENDIX  B 


Characteristic  Boundary  Conditions;  P  &  U  Specified  (McKenna) 
The  characteristic  variables  are 


qx  -  p  -  p/c; 


Q2  =  P  +  Pocou 


Q3  -  P  -  PoCoU 


In  terms  of  primitive  variables  they  are 


q2-q3 


q2-q3 


P  -  + 


2p0C0 

Q,  +  Q-: 

1  +  2CT 


These  values  were  derived  by  solving  the  three  equations  (above)  for 
three  unknowns. 

This  information  is  now  applied  at  the  subsonic  (exit)  boundary. 
Subsonic  Outflow: 


The  boundary  conditions  are  mathematically  represented  as 


constant  =  K_ 


Refer  now  to  the  primitive  variables  written  in  terms  of  the 
characteristics.  The  velocity,  for  example  is  written  as 


^2  ~  ^3 

U  *  -4 - - 

2p(jC0 

This  term  contains  two  characteristics  -  the  second  and  third.  McKenna 
reasoned  that  is  constant  while  and  (at  the  exit)  are  allowed 
to  float.  Therefore  in  order  to  set  the  exit  boundary  condition,  Qj  and 
should  be  expressed  in  terms  of  the  primitive  variables.  is  held 
constant  at  K^.  The  velocity  term  would  then  become 

^2  Q3 

'  P  +  P0C0U  -  k3 

The  other  terms  are 


P  - 


(P  ♦  PqCqU)  +  K3 


The  variables  are  numerically  implemented  into  the  program  as: 

K- 


Pn  »  £a£a. 

j  2 


*  Uj-1  *  ^7] 

Pj  -  pj-i  *  [0^7  -  pfo  *  °j-ij 
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^  - 


UU  V  -  A 


Where  pfl  &  C0  are  evaluated  at  the  j  point  for  the  n  -  1  time-step 
Kj  is  a  constant  represented  as 


Pm  - 


CnPnUM 
o  o  00 


Here  and  are  the  values  at  the  exit. 


APPENDIX  C 


Characteristic  Boundary  Condition:  P  specified  (Steger) 
Subsonic  Outflow: 

Steger  represented  the  outflow  boundary  as 


P  -  Pea 


This  can  be  rewritten  as 


Qi.  ■  qi.  , 
J  J-l 


Q2.  "  Q2.  , 

J  J-l 


P  -  Pa 


In  terms  of  the  primitive  variables  they  are  written  as 


p.  -  P./C*  -  p.  -  -  P.  ,/Cn 

j  j  0  j-l  j-l  0 


P.  +  p„C„U.  «  P.  .  +  pnCn  U.  . 
1  K0  0  1  1-1  0  1-1 


Now  solving  for  the  density  and  velocity  we  get 


Pj  "  pj-i  +  cj  ^Pco  '  Pj-l) 

"j  -  Vi ♦  ok  (Pi-i  •  P~> 


pj  =  p“ 


Numerically  these  equations  are  represented  as 


n  n  1  /n  \ 

Pj  =  Pj-l  +  c7  (P°°  ”  Pj-i) 


„n  ..n  1  .  __n 

U.  -  U.  ,  +  (P.  , 

J  J"1  Poco  11  J-1 


-  Pa 


P"  B  P„ 


p„  and  C.  are  the  same  as  before. 
Ko  o 
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APPENDIX  D 


Sample  Calculations  for  the  Initial  Shock  Location: 

P0  =  10,000  lbf/ft2 
Given : 


p0  =  .00883  lbfsec2ft^ 


The  values  of  various  points  are: 
PT  02 

(X  =  .4)  A2  =  1.0514 


M2  =  1.5 


.4“  1*176  J-  =  *2724  £■  =  .3950 


A* 


A*  -  .8941  P2  *  2724 


U  =  1568.45 
2 


Po 

P, 


-  C034  C 


PT  01 


(X  =  0.0) 


1  1.0512 

A*  "  .894) 


=  1.1757  M  =  1.4996 


.2726  =  .3951  C  =  1045.9 


Ux  =  1568.38  P  -  2726  pj  =  .0034 


-  1045. 
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The  shock  is  at  point  16 
PT  #16 


(X  =  6.0) 


A16  1.6284 

A*  "  .8941 


1.821  Mlfi 


2.09 


f  -  .1111 


=  .2081  C  =  920.04 
po 


U16  =  1922‘9  P16  =  1111  Pl6  =  -0' 

After  the  shock 


M  -  .5628 

y 

p  2  _ 
pl 

4.929 

=  2.798 

Ci 

P0V  _ 
p°x 

.6789 

P 1  6  = 
P02 

.1636 

U  =  687.3 

y 

P  = 

y 

5476.12 

p  =  .0051 
y 

p 

P  0  y  =  ‘ 

.8065 

.  P  —  ^ 

POy 

8576 

P0  »  6789 

y 

POy  1 

=  .0059 

A  =  1.2363 

y 

M2  s 


Ae  _  A2 5  A*  Ai  b  Ay  _  A2 5 
Ae  A*  Ai  6  Ay  A*^  A* 

AE  1.7445  .1.1  1.2363 

A*  =  .8941  1.821  1  '  1 

=  .508  ~~  =  .8385  =  .8817  CZ5 

*  0y  POy 


1.3248 

=  1228.6 


PT  #25 

(X  *=  9.6)  U25  =  624  P25  -  5692.6  p25  =  .00: 

PT  #26 


(X  =  10.0) 


& 


1.7447 
.8941 


1 

1.821 


1.2363 

1 


1.3249 


(same  as  PT  if 25) 


u 

P 

£ 

M 

1568.38 

2726 

.0034 

1.4996 

1568.45 

2724 

.0034 

1.5 

1922.9 

1111 

.0018 

2.09 

687.3 

5476.12 

.0051 

.5628 

623.99 

5692.6 

.0052 

.508 

623.0 

5692.9 

.0052 

.5080 

APPENDIX  E 


Theoretical  Values  for  Moving  the  Shock  Upstream  and  Downstream 

Upstream;  (From  the  original  shock)  Choose  point  #12  for  the 
new  shock  location 

PT  #12  =  =  1,3904  Mi2  =  1.7540 

~  =  .1867  =  .3016 

P  o 

t  Mt 

Pi  2  =  1867  =  .0027  U  -  1725.8 

ML.  =  .6271  tt*-  =  3.422  =  .8328  ^-=2.28 

"  Pl2  Pt12  P I 2 

P  =  6388.9  p  =  .0063  U  =  753.2  P_  =  8328  p_  =  .0' 
y  y  y  ty  ^ty 

The  exit  conditions  are  now  calculated 


Downstream:  (From  the  original  shock)  Choose  point  it 20 
new  downstream  shock 

PT  #20  =  1.9398  M20  =  2.1827  My  = 

=  .1919  =  .6408  p20  =  .0017  |^-  = 

fci  1 

M  =  .5521  4"  =  1*2519  =  *8625  P_  =  6408 

A*  pt  ty 


for  the 

.5521 

.8997 


The  exit  conditions  are  now  calculated 


APPENDIX  F  -  Program  Listing 
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20.  ABSTRACT  ( Continue  on  reverse  side  If  necessary  and  Identify  by  block  number) 

There  are  basically  two  techniques  used  to  solve  the  Navier-Stokes 
equations  for  fluid  flow.  These  techniques  are  implicit  and  explicit  methods 
Both  methods  along  with  characteristic  boundary  conditions  for  solving  quasi- 
one-dimensional  nozzle  flow  are  presented.  Two  types  of  characteristic 
boundary  conditions,  those  of  Bteger  and  McKenna,  were  tested  for  each 
scheme.  Solutions  of  isentropic  supersonic  flow  and  flow  with  shocks  were 
obtained  for  a  diverging  nozzle.  'The  behavior  of  each  boundary  condition  on 
both  implicit  and  explicit  schemes  were  the  sane.  They  deviated  from  the 
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theoretical  values  by  less  than  one  percent.  A  rest  to  determine  the 
utility  of  each  scheme  was  run  by  a} lowing  the  exit  boundary  conditions 
to  change  in  the  hope  of  forcing  the  shock  to  move  upstream  or  down¬ 
stream.  The  shock  would  not  move  in  the  implicit  scheme  for  eitner 
boundary  condition.  The  explicit  scheme  could  move  the  shock,  but  only 
when  Steger's  boundary  condition  was  used.  (This  is  the  one  which 
specified  only  pressure.) 


UI.TCLAS3IFISD 


SECURITY  CLASSIFICATION  OF  THIS  PAGEfHTi.n  D.r.  Enf.r.rf) 


